
%% initialize
clear all;
close all;

rng(328120)

Nsim = 500;

loadstarting = 0;
loadestimates = 0;

%% read in data

load '../../01_data/05_temp/mixed_logit_rc6_totVA_DISAD_Ndis';

% convert 20% of applications to non-applications based on utility index
[Xdraws_best] = func_likely_draws_current(IJitidx, IJchoiceset,IJcurrentFirst, Xdraws,Xmat,thetahat,y_logistic);
    
    % redo Xdraws_best to be as long as everything else
    [Ci, IAi, ICi]=unique(IJitidx(IJchoiceset==1));
    Xdraws_unique=Xdraws_best(IAi,:);
    Xdraws_best2=Xdraws_unique(IJitidx,:);
    
    % find predicted values in each match (teacher preferences)
    
    IJipref=[Xmat Xdraws_best2]*thetahat;
    
epsilon = evrnd(0,1,size(IJipref,1),1);   

IJiprefeps = IJipref+epsilon;

IJpercentiles=cell2mat(accumarray(IJitidx,IJiprefeps,[], @(x) {prctile(x,[20])}));

IJpercentileslong = IJpercentiles(IJitidx);

top80 = IJiprefeps>=IJpercentileslong;

numapp = accumarray(IJitidx,IJapplied,[],@sum);
numtop80 = accumarray(IJitidx,top80.*IJapplied,[],@sum);

numapplong = numapp(IJitidx);
numtop80long = numtop80(IJitidx);

IJappliedorig = IJapplied;

IJapplied(top80==0 & numtop80long>0) = 0;

f =  @(x)func_mixed_logit_current(x,IJapplied(IJchoiceset==1),Xmat(IJchoiceset==1,:),...
    Xdraws(IJchoiceset==1,:,:),IJitidx(IJchoiceset==1),Xmat(IJcurrentFirst==1,:),...
    Xdraws(IJcurrentFirst==1,:,:),IJitidx(IJcurrentFirst==1),0);


tic
thetahat = fminunc(f,thetastart,optionsgrad);
toc

 save -v7.3 '../../01_data/05_temp/mixed_logit_rc6_totVA_convert'
 
